


###########################################################################
###########################################################################
## ZIMBABWE ANALYSIS - VIGNETTE EXPERIMENT
## 2019 Mar 25
## Lauren Young
###########################################################################
###########################################################################


rm(list=ls())

library(data.table)
library(stargazer)
library(foreign)
library(readstata13)
library(gtools)
library(ggplot2)
library(reshape2)
library(interplot)
library(psych)
library(xtable)

set.seed(11010)

setwd('/Users/leyou/Box/Zimbabwe/self_efficacy/01_RepFiles/')

source('02_Code/clusterfxn.R')

st <- function(x) {
  (x - mean(x,na.rm=T))/sd(x,na.rm=T)
}

se <- function(x){sd(x, na.rm = T)/sqrt(length(na.omit(x)))}



#####
## read in data
#####

## prevalence of state-sponsored election violence
pts <- read.dta13("01_Data/01_StrausPevehouseTaylor_AEVD_2017/PerilsPluralism1.dta")


## afrobarometer data
afro <- read.csv('01_Data/04_Afrobarometer/afro_clean_r5.csv')


## read in 2015 survey data
dat <- read.csv('01_Data/02_2015Survey/2015_indiv.csv')
con <- read.csv('01_Data/02_2015Survey/2015_stacked.csv')



#####
## Table 1: SEF, anger, and dissent 
#####

spec0 <- "con.con_wa_rec + con.con_wb_sim + con.con_wc_same + con.con_wd_sev + con.con_we_cred + se_pc_st + treat_assign + community + surveyor_id"
spec1 <- " + party_close_st + female + age_st + education_st + assets_pc_urb_st + assets_pc_rur_st "
spec2 <- "+ activism_pc_st + trauma_pc_st"
spec3 <- " + con.con_wa_rec:se_pc_st + con.con_wb_sim:se_pc_st + con.con_wc_same:se_pc_st + con.con_wd_sev:se_pc_st + con.con_we_cred:se_pc_st"

## probability of rally

mod4a <- lm(as.formula(paste0('act_st ~ ', spec0)), data = con)
mod4a$se <- coeftest.cluster(con, mod4a, 'id')[,2]
mod4b <- lm(as.formula(paste0('act_st ~ ', spec0, spec3, spec1)), data = con)
mod4b$se <- coeftest.cluster(con, mod4b, 'id')[,2]
mod4c <- lm(as.formula(paste0('act_st ~ ', spec0, spec3, spec1, spec2)), data = con)
mod4c$se <- coeftest.cluster(con, mod4c, 'id')[,2]
mod4d <- lm(as.formula(paste0('act_st ~ ', spec0, spec1, spec2, spec3)), data = con)
mod4d$se <- coeftest.cluster(con, mod4d, 'id')[,2]


## anger - fear
mod2a <- lm(as.formula(paste0('ratio_af_st ~ ', spec0)), data = con)
mod2a$se <- coeftest.cluster(con, mod2a, 'id')[,2]
mod2b <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, spec3, spec1)), data = con)
mod2b$se <- coeftest.cluster(con, mod2b, 'id')[,2]
mod2c <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, spec3, spec1, spec2)), data = con)
mod2c$se <- coeftest.cluster(con, mod2c, 'id')[,2]
mod2d <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, spec1, spec2, spec3)), data = con)
mod2d$se <- coeftest.cluster(con, mod2d, 'id')[,2]

stargazer(mod4a, mod4b, mod4c, mod2a, mod2b, mod2c, 
          se = list(mod4a$se, mod4b$se, mod4c$se, mod2a$se, mod2b$se, mod2c$se),
          covariate.labels = c('Time to Election', 'Non-Activist Victim', 'Same Province',
                               'Violence Severity', 'Source Credibility',
                               'Self Efficacy', 'Closeness to Party',
                               'Female', 'Age',
                               'Education', 'Urban Assets', 'Rural Assets',
                               'Past Activism', 'Past Trauma',
                               'Self Efficacy X Time to Election',
                               'Self Efficacy X Non-Activism Victim',
                               'Self Efficacy X Same Province',
                               'Self Efficacy X Violence Severity',
                               'Self Efficacy X Source Credibility',
                               'Constant'),
          omit = c('treat_assign', 'community', 'surveyor_id'),
          keep.stat = c('n', 'rsq'),
          no.space = T, digits = 2, star.cutoffs = c(.05, .01))




#####
## Figure 1: Plot marginal effects
#####


## self efficacy X severity - dissent

p <- interplot(mod4c, var1 = 'con.con_wd_sev', var2 = 'se_pc_st', hist = T)

pdf('03_Figures/dissent_marg_sef.pdf')
p + 
  xlab('Self-Efficacy') +
  ylab('Marginal Effect of Violence Severity') +
  theme_minimal() +
  theme(text = element_text(size=18), 
        # axis.text.x = element_text(angle=45, vjust=0.95, hjust=1),
        legend.position = "bottom",
        panel.spacing = unit(0.5, "lines")) 
dev.off()


## self efficacy X severity - ratio anger-fear

p <- interplot(mod2c, var1 = 'con.con_wd_sev', var2 = 'se_pc_st', hist = T)

pdf('03_Figures/ratio_marg_sef.pdf')
p + 
  xlab('Self-Efficacy') +
  ylab('Marginal Effect of Violence Severity') +
  theme_minimal() +
  theme(text = element_text(size=18), 
        # axis.text.x = element_text(angle=45, vjust=0.95, hjust=1),
        legend.position = "bottom",
        panel.spacing = unit(0.5, "lines")) 
dev.off()




#####
## Appendix table D.1: sample characteristics
#####


## subset AB and WB to only opposition supporters
afro_all <- afro
afro <- afro[afro$party_zanu==0 & afro$party_none==0,]
## NB: code to replicate the world bank analyses is not made available for replication because I do not have permission to share the dataset.
# wb_all <- wb
# wb <- wb[wb$party_which!='ZANU' & wb$party_any==1,]


tab <- data.frame('label' = c('age', 'education', 'female', 'food', 'income', 'urban',
                              'party_close', 'sef_means', 'sef_coping', 'sef_resource', 'sef_prob'),
                  'mean_ly' = NA, 'sd_ly' = NA, 'n_ly' = NA,
                  'mean_ab' = NA, 'sd_ab' = NA, 'n_ab' = NA,
                  'mean_wb' = NA, 'sd_wb' = NA, 'n_wb' = NA)

afro$econinc <- afro$nocash+1
afro$econfood <- afro$hunger+1

tab[1, c(2:10)] <- c(mean(dat$age, na.rm=T), sd(dat$age,na.rm=T), length(na.omit(dat$age)),
                     mean(afro$age, na.rm=T), sd(afro$age,na.rm=T), length(na.omit(afro$age)),
                     # mean(wb$age, na.rm=T), sd(wb$age,na.rm=T), length(na.omit(wb$age))
                     NA, NA, NA)
tab[2, c(2:10)] <- c(mean(dat$education, na.rm=T), sd(dat$education,na.rm=T), length(na.omit(dat$education)),
                     mean(afro$educ, na.rm=T), sd(afro$educ,na.rm=T), length(na.omit(afro$educ)),
                     # mean(wb$education, na.rm=T), sd(wb$education,na.rm=T), length(na.omit(wb$education))
                     NA, NA, NA)
tab[3, c(2:10)] <- c(mean(dat$female, na.rm=T), sd(dat$female,na.rm=T), length(na.omit(dat$female)),
                     mean(afro$female, na.rm=T), sd(afro$female,na.rm=T), length(na.omit(afro$female)),
                     NA, NA, NA)
                     # mean(wb$female, na.rm=T), sd(wb$female,na.rm=T), length(na.omit(wb$female)))
tab[4, c(2:10)] <- c(mean(dat$econfood, na.rm=T), sd(dat$econfood,na.rm=T), length(na.omit(dat$econfood)),
                     mean(afro$econfood, na.rm=T), sd(afro$econfood,na.rm=T),  length(na.omit(afro$econfood)),
                     NA, NA, NA)
                     # mean(wb$econfood, na.rm=T), sd(wb$econfood,na.rm=T), length(na.omit(wb$econfood)))
tab[5, c(2:10)] <- c(mean(dat$econinc, na.rm=T), sd(dat$econinc,na.rm=T), length(na.omit(dat$econinc)),
                     mean(afro$econinc, na.rm=T), sd(afro$econinc,na.rm=T), length(na.omit(afro$econinc)),
                     NA, NA, NA)
# mean(wb$econinc, na.rm=T), sd(wb$econinc,na.rm=T), length(na.omit(wb$econinc)))
tab[6, c(2:10)] <- c(mean(dat$urban, na.rm=T), sd(dat$urban,na.rm=T), length(na.omit(dat$urban)),
                     mean(afro$urb, na.rm=T), sd(afro$urb,na.rm=T), length(na.omit(afro$urb)),
                     NA, NA, NA)
tab[7, c(2:10)] <- c(mean(dat$party_close, na.rm=T), sd(dat$party_close,na.rm=T), length(na.omit(dat$party_close)),
                     # mean(afro$party_close, na.rm=T), sd(afro$party_close,na.rm=T), length(na.omit(afro$party_close)),
                     NA, NA, NA,
                     # mean(wb$party_close, na.rm=T), sd(wb$party_close,na.rm=T), length(na.omit(wb$party_close)))
                     NA, NA, NA)
tab[8, c(2:10)] <- c(mean(dat$sef_means, na.rm=T), sd(dat$sef_means,na.rm=T), length(na.omit(dat$sef_means)),
                     # mean(afro$party_close, na.rm=T), sd(afro$party_close,na.rm=T), length(na.omit(afro$party_close)),
                     NA, NA, NA,
                     # mean(wb$sef_means, na.rm=T), sd(wb$sef_means,na.rm=T), length(na.omit(wb$sef_means)))
                     NA, NA, NA)
tab[9, c(2:10)] <- c(mean(dat$sef_coping, na.rm=T), sd(dat$sef_coping,na.rm=T), length(na.omit(dat$sef_coping)),
                     # mean(afro$party_close, na.rm=T), sd(afro$party_close,na.rm=T), length(na.omit(afro$party_close)),
                     NA, NA, NA,
                     # mean(wb$sef_coping, na.rm=T), sd(wb$sef_coping,na.rm=T), length(na.omit(wb$sef_coping)))
                     NA, NA, NA)
tab[10, c(2:10)] <- c(mean(dat$sef_resource, na.rm=T), sd(dat$sef_resource,na.rm=T), length(na.omit(dat$sef_resource)),
                      # mean(afro$party_close, na.rm=T), sd(afro$party_close,na.rm=T), length(na.omit(afro$party_close)),
                      NA, NA, NA,
                      # mean(wb$sef_resource, na.rm=T), sd(wb$sef_resource,na.rm=T), length(na.omit(wb$sef_resource)))
                      NA, NA, NA)
tab[11, c(2:10)] <- c(mean(dat$sef_prob, na.rm=T), sd(dat$sef_prob,na.rm=T), length(na.omit(dat$sef_prob)),
                      # mean(afro$party_close, na.rm=T), sd(afro$party_close,na.rm=T), length(na.omit(afro$party_close)),
                      NA, NA, NA,
                      # mean(wb$sef_prob, na.rm=T), sd(wb$sef_prob,na.rm=T), length(na.omit(wb$sef_prob)))
                      NA, NA, NA)

rownames(tab) <- c('Age', 'Education', 'Female', 'Subjective Poverty: Food', 
                   'Subjective Poverty: Income', 'Urban', 'Closeness to Party', 
                   'SEF: Means', "SEF: Coping", 'SEF: Resources', 'SEF: Problems')
tab$label=NULL

xtable(tab, digits = 2)




#####
## Appendix Figure E.1: self efficacy validation
#####

vars <- c('se.se1', 'se.se2', 'se.se3', 'se.se4', 'se.se5',
          'se.se6', 'se.se7', 'se.se8', 'se.se9', 'se.se10')
psych::alpha(dat[,vars])

cormat <- round(cor(dat[,vars], use = 'complete.obs'),2)
cormat[lower.tri(cormat)] <- NA


cormat <- na.omit(melt(cormat, na.rm=T))
cormat$Var1 <- gsub("se.se", "SE ", as.character(cormat$Var1))
cormat$Var2 <- gsub("se.se", "SE ", as.character(cormat$Var2))
cormat$Var1 <- factor(cormat$Var1, levels = paste0('SE ', seq(1,10,1)))
cormat$Var2 <- factor(cormat$Var2, levels = paste0('SE ', seq(1,10,1)))


pdf("03_Figures/sef_heatmap.pdf")
ggplot(data = cormat, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile(color = "white") +
  geom_text(aes(label = value), color = "black", size = 4) +
  theme_minimal() +
  theme(text = element_text(size=14)) +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson\nCorrelation") +
  coord_fixed()
dev.off()



#####
## appendix table E.1: validation of scenario reactions
#####

spec1 <- '+ surveyor_id + community + treat_assign'
spec2 <- "+ party_close_st + education_st + female + age_st + assets_pc_urb_st + assets_pc_rur_st"


mod1a <- lm(as.formula(paste0('act_st ~ ', 'prob_act_st', spec1)), data = con)
mod1a$se <- coeftest.cluster(con, mod1a, 'id')[,2]
mod1b <- lm(as.formula(paste0('act_st ~ ', 'prob_act_st', spec1, spec2)), data = con)
mod1b$se <- coeftest.cluster(con, mod1b, 'id')[,2]

mod2a <- lm(as.formula(paste0('act_st ~ ', 'wristband', spec1)), data = con[con$wristband_real==1,])
mod2a$se <- coeftest.cluster(con[con$wristband_real==1,], mod2a, 'id')[,2]
mod2b <- lm(as.formula(paste0('act_st ~ ', 'wristband', spec1, spec2)), data = con[con$wristband_real==1,])
mod2b$se <- coeftest.cluster(con[con$wristband_real==1,], mod2b, 'id')[,2]

mod3a <- lm(as.formula(paste0('act_st ~ ', 'activism_pc_st', spec1)), data = con)
mod3a$se <- coeftest.cluster(con, mod3a, 'id')[,2]
mod3b <- lm(as.formula(paste0('act_st ~ ', 'activism_pc_st', spec1, spec2)), data = con)
mod3b$se <- coeftest.cluster(con, mod3b, 'id')[,2]

stargazer(mod1a, mod1b, mod3a, mod3b, mod2a, mod2b,
          covariate.labels = c('Opposition Propensity (Hypothetical)',
                               'Past Activism (Self Reported)',
                               'Wristband (Behavioral)',
                               'Closeness to Party',
                               'Female', 'Age', 'Education',
                               'Urban Assets', 'Rural Assets', 'Constant'),
          se = list(mod1a$se, mod1b$se, mod3a$se, mod3b$se, mod2a$se, mod2b$se),
          omit = c('treat_assign', 'community', 'surveyor_'), no.space = T, digits = 2,
          star.cutoffs = c(.05, .01))



#####
## Appendix Table F.1: demog correlates of self efficacy
#####


spec1 <- "female + age_st + education_st + assets_pc_urb_st + assets_pc_rur_st + community"
spec2 <- "+ party_close_st"
spec2wb <- "+ party_close_st + socgrp_women + socgrp_union + socgrp_savings + socgrp_farming + socgrp_relig + socgrp_residents + socgrp_other"
spec3 <- "+ activism_pc_st"
spec4 <- "+ trauma_pc_st"

mod1a <- lm(as.formula(paste0('se_pc_st ~ ', spec1, "+ treat_assign")), data = dat)
mod1b <- lm(as.formula(paste0('se_pc_st ~ ', spec1, "+ treat_assign", spec2)), data = dat)
mod1c <- lm(as.formula(paste0('se_pc_st ~ ', spec1, "+ treat_assign", spec2, spec3)), data = dat)
mod1d <- lm(as.formula(paste0('se_pc_st ~ ', spec1, "+ treat_assign", spec2, spec3, spec4)), data = dat)

# mod2a <- lm(as.formula(paste0('se_pc_st ~ ', spec1)), data = wb)
# mod2b <- lm(as.formula(paste0('se_pc_st ~ ', spec1, spec2wb)), data = wb)
# mod2c <- lm(as.formula(paste0('se_pc_st ~ ', spec1, spec2wb, spec3)), data = wb)
# mod2d <- lm(as.formula(paste0('se_pc_st ~ ', spec1, spec2wb, spec3, spec4)), data = wb)


stargazer(mod1a, mod1b, mod1c, mod1d, #mod2a, mod2b, mod2c, mod2d, 
          covariate.labels = c('Female', 'Age', 'Education',
                               'Urban Assets', 'Rural Assets',
                               'Closeness to Party',
                               'Activism Experience', 'Violence Exposure',
                               # 'Membership: Womens Group', 'Membership: Union',
                               # 'Membership: Savings Group', 'Membership: Farming Group',
                               # 'Membership: Religious Group', 'Membership: Residents Assoc.',
                               # 'Membership: Other Group',
                               'Constant'),
          omit = c('treat_assign', 'community'),
          keep.stat = c('n', 'rsq'),
          no.space = T, digits = 2,
          star.cutoffs = c(.05, .01))




#####
## Appendix Table G.1: Disaggregate anger and fear
#####

spec0 <- "con.con_wa_rec + con.con_wb_sim + con.con_wc_same + con.con_wd_sev + con.con_we_cred + se_pc_st + treat_assign + community + surveyor_id"
spec1 <- " + party_close_st + female + age_st + education_st + assets_pc_urb_st + assets_pc_rur_st "
spec2 <- "+ activism_pc_st + trauma_pc_st"
spec3 <- " + con.con_wa_rec:se_pc_st + con.con_wb_sim:se_pc_st + con.con_wc_same:se_pc_st + con.con_wd_sev:se_pc_st + con.con_we_cred:se_pc_st"

## anger

mod4a <- lm(as.formula(paste0('anger_st ~ ', spec0)), data = con)
mod4a$se <- coeftest.cluster(con, mod4a, 'id')[,2]
mod4b <- lm(as.formula(paste0('anger_st ~ ', spec0, spec3, spec1)), data = con)
mod4b$se <- coeftest.cluster(con, mod4b, 'id')[,2]
mod4c <- lm(as.formula(paste0('anger_st ~ ', spec0, spec3, spec1, spec2)), data = con)
mod4c$se <- coeftest.cluster(con, mod4c, 'id')[,2]


## fear

mod2a <- lm(as.formula(paste0('fear_st ~ ', spec0)), data = con)
mod2a$se <- coeftest.cluster(con, mod2a, 'id')[,2]
mod2b <- lm(as.formula(paste0('fear_st ~ ', spec0, spec3, spec1)), data = con)
mod2b$se <- coeftest.cluster(con, mod2b, 'id')[,2]
mod2c <- lm(as.formula(paste0('fear_st ~ ', spec0, spec3, spec1, spec2)), data = con)
mod2c$se <- coeftest.cluster(con, mod2c, 'id')[,2]

## Table G.1

stargazer(mod4a, mod4b, mod4c, mod2a, mod2b, mod2c, 
          se = list(mod4a$se, mod4b$se, mod4c$se, mod2a$se, mod2b$se, mod2c$se),
          covariate.labels = c('Time to Election', 'Non-Activist Victim', 'Same Province',
                               'Violence Severity', 'Source Credibility',
                               'Self Efficacy', 'Closeness to Party',
                               'Female', 'Age',
                               'Education', 'Urban Assets', 'Rural Assets',
                               'Past Activism', 'Past Trauma',
                               'Self Efficacy X Time to Election',
                               'Self Efficacy X Non-Activism Victim',
                               'Self Efficacy X Same Province',
                               'Self Efficacy X Violence Severity',
                               'Self Efficacy X Source Credibility',
                               'Constant'),
          omit = c('treat_assign', 'community', 'surveyor_id'),
          keep.stat = c('n', 'rsq'),
          no.space = T, digits = 2, star.cutoffs = c(.05, .01))



#####
## Appendix G.2: test of overreporting
#####

## table G.2

tab <- table(con$act_st, con$wristband, useNA = 'ifany')
tab2 <- prop.table(tab, margin = 2)
rownames(tab) <- c('Not at all', 'A little bit', 'Somewhat', 'Very', 'Sure', 'Missing')


tab <- data.frame('Refused' = tab[,1], 'Proportion Ref' = tab2[,1],
                  'Accepted' = tab[,2], 'Proportion Acc' = tab2[,2],
                  'Missing' = tab[,3], 'Proportion Miss' = tab2[,3])
xtable(tab)


## table G.3

con$hyp_only0 <- ifelse(con$wristband==0 & con$act_st>0, 1, 0)
con$hyp_only0[is.na(con$wristband)==T] <- NA
con$hyp_only1 <- ifelse(con$wristband==0 & con$act_st>.1, 1, 0)
con$hyp_only1[is.na(con$wristband)==T] <- NA
con$hyp_only2 <- ifelse(con$wristband==0 & con$act_st>1, 1, 0)
con$hyp_only2[is.na(con$wristband)==T] <- NA

spec1 <- 'se_pc_st + treat_assign + community + surveyor_id'
spec2 <- "+ party_close_st + education_st + female + age_st + assets_pc_urb_st + assets_pc_rur_st"
spec3 <- "+ activism_pc_st + trauma_pc_st"

mod1a <- lm(as.formula(paste0('hyp_only0 ~ ', spec1)), data = con)
mod1a$se <- coeftest.cluster(con, mod1a, 'id')[,2]
mod1b <- lm(as.formula(paste0('hyp_only0 ~ ', spec1, spec2)), data = con)
mod1b$se <- coeftest.cluster(con, mod1b, 'id')[,2]
mod1c <- lm(as.formula(paste0('hyp_only0 ~ ', spec1, spec2, spec3)), data = con)
mod1c$se <- coeftest.cluster(con, mod1c, 'id')[,2]

mod2a <- lm(as.formula(paste0('hyp_only1 ~ ', spec1)), data = con)
mod2a$se <- coeftest.cluster(con, mod2a, 'id')[,2]
mod2b <- lm(as.formula(paste0('hyp_only1 ~ ', spec1, spec2)), data = con)
mod2b$se <- coeftest.cluster(con, mod2b, 'id')[,2]
mod2c <- lm(as.formula(paste0('hyp_only1 ~ ', spec1, spec2, spec3)), data = con)
mod2c$se <- coeftest.cluster(con, mod2c, 'id')[,2]

mod3a <- lm(as.formula(paste0('hyp_only2 ~ ', spec1)), data = con)
mod3a$se <- coeftest.cluster(con, mod3a, 'id')[,2]
mod3b <- lm(as.formula(paste0('hyp_only2 ~ ', spec1, spec2)), data = con)
mod3b$se <- coeftest.cluster(con, mod3b, 'id')[,2]
mod3c <- lm(as.formula(paste0('hyp_only2 ~ ', spec1, spec2, spec3)), data = con)
mod3c$se <- coeftest.cluster(con, mod3c, 'id')[,2]

stargazer(mod1a, mod1b, mod1c, mod2a, mod2b, mod2c, mod3a, mod3b, mod3c,  
          covariate.labels = c('Self-Efficacy',
                               'Closeness to Party', 'Education', 'Female',
                               'Age', 'Urban Assets', 'Rural Assets',
                               'Activism Experience', 'Violence Exposure',
                               'Constant'),
          se = list(mod1a$se, mod1b$se, mod1c$se, mod2a$se, mod2b$se, mod2c$se, mod3a$se, mod3b$se, mod3c$se),
          omit = c('treat_assign', 'community', 'surveyor_'), no.space = T, digits = 2,
          keep.stat=c('n', 'rsq'), star.cutoffs = c(.05, .01))




#####
## Appendix Table G.3: SEF, anger, and dissent - disagg by party
#####

spec0 <- "con.con_wa_rec + con.con_wb_sim + con.con_wc_same + con.con_wd_sev + con.con_we_cred + se_pc_st + treat_assign + community + surveyor_id"
spec1 <- " + party_close_st + female + age_st + education_st + assets_pc_urb_st + assets_pc_rur_st "
spec2 <- "+ activism_pc_st + trauma_pc_st"
spec3 <- " + con.con_wa_rec:se_pc_st + con.con_wb_sim:se_pc_st + con.con_wc_same:se_pc_st + con.con_wd_sev:se_pc_st + con.con_we_cred:se_pc_st"


## probability of rally

mod4a <- lm(as.formula(paste0('act_st ~ ', spec0, "+ party_which_mdc + se_pc_st:party_which_mdc")), data = con)
mod4a$se <- coeftest.cluster(con, mod4a, 'id')[,2]
mod4b <- lm(as.formula(paste0('act_st ~ ', spec0, "+ party_which_mdc + se_pc_st:party_which_mdc", spec3, spec1)), data = con)
mod4b$se <- coeftest.cluster(con, mod4b, 'id')[,2]
mod4c <- lm(as.formula(paste0('act_st ~ ', spec0, "+ party_which_mdc + se_pc_st:party_which_mdc", spec3, spec1, spec2)), data = con)
mod4c$se <- coeftest.cluster(con, mod4c, 'id')[,2]


## anger - fear

mod2a <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, "+ party_which_mdc + se_pc_st:party_which_mdc")), data = con)
mod2a$se <- coeftest.cluster(con, mod2a, 'id')[,2]
mod2b <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, "+ party_which_mdc + se_pc_st:party_which_mdc", spec3, spec1)), data = con)
mod2b$se <- coeftest.cluster(con, mod2b, 'id')[,2]
mod2c <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, "+ party_which_mdc + se_pc_st:party_which_mdc", spec3, spec1, spec2)), data = con)
mod2c$se <- coeftest.cluster(con, mod2c, 'id')[,2]


## table
stargazer(mod4a, mod4b, mod4c, mod2a, mod2b, mod2c, 
          se = list(mod4a$se, mod4b$se, mod4c$se, mod2a$se, mod2b$se, mod2c$se),
          covariate.labels = c('Time to Election', 'Non-Activist Victim', 'Same Province',
                               'Violence Severity', 'Source Credibility',
                               'Self Efficacy', 'MDC Supporter', 'Closeness to Party',
                               'Female', 'Age',
                               'Education', 'Urban Assets', 'Rural Assets',
                               'Past Activism', 'Past Trauma',
                               'Self Efficacy X MDC Supporter',
                               'Self Efficacy X Time to Election',
                               'Self Efficacy X Non-Activism Victim',
                               'Self Efficacy X Same Province',
                               'Self Efficacy X Violence Severity',
                               'Self Efficacy X Source Credibility',
                               'Constant'),
          omit = c('treat_assign', 'community', 'surveyor_id'),
          keep.stat = c('n', 'rsq'),
          no.space = T, digits = 2, star.cutoffs = c(.05, .01))




#####
## Table (not reported): Categorical scenario treatments
#####

con$con.con_wa <- factor(con$con.con_wa, levels = c('year', 'month', 'day'))
con$con.con_wb <- factor(con$con.con_wb, levels = c('parliamentary candidate', 'council candidate', 'organizer', 'voter', 'friend'))
con$con.con_wd <- factor(con$con.con_wd, levels = c('threatened', 'beaten', 'abducted', 'killed'))
con$con.con_we <- factor(con$con.con_we, levels = c('a ZANU-PF activist', 'an opposition activist', 'your friend'))

spec0 <- "con.con_wa + con.con_wb + con.con_wc_same + con.con_wd + con.con_we + treat_assign + community + surveyor_id"


mod1a <- lm(as.formula(paste0('ratio_af_st ~ ', spec0)), data = con)
mod1a$se <- coeftest.cluster(con, mod1a, 'id')[,2]

mod2a <- lm(as.formula(paste0('fear_st ~ ', spec0)), data = con)
mod2a$se <- coeftest.cluster(con, mod2a, 'id')[,2]

mod3a <- lm(as.formula(paste0('anger_st ~ ', spec0)), data = con)
mod3a$se <- coeftest.cluster(con, mod3a, 'id')[,2]

mod4a <- lm(as.formula(paste0('act_st ~ ', spec0)), data = con)
mod4a$se <- coeftest.cluster(con, mod4a, 'id')[,2]


#####
## Figure G.1: Categorical scenario treatments
#####

## anger and fear and rally

p <- data.frame('Type' = c(rep('time', 9), rep('activism', 15), rep('loc', 6),
                           rep('sev', 12), rep('cred', 9)),
                'Outcome' = c(rep(c("Fear", "Anger", "Rally"), 17)),
                'Characteristic' = c(rep('Time to Election: Year', 3),
                                     rep('Time to Election: Month', 3), 
                                     rep('Time to Election: Day', 3), 
                                     rep('Victim Activism: Parliamentary candidate', 3),
                                     rep('Victim Activism: Council candidate', 3), 
                                     rep('Victim Activism: Local mobilizer', 3), 
                                     rep('Victim Activism: Voter', 3),
                                     rep('Victim Activism: Friend', 3), 
                                     rep('Same Province: No', 3), rep('Same Province: Yes', 3), 
                                     rep('Violence Severity: Threatened', 3), 
                                     rep('Violence Severity: Beaten', 3),
                                     rep('Violence Severity: Abducted', 3), 
                                     rep('Violence Severity: Killed', 3), 
                                     rep('Source Credibility: Ruling party activist', 3),
                                     rep('Source Credibility: Opposition activist', 3),
                                     rep('Source Credibility: Your friend', 3)),
                'diffs' = c(0, 0, 0,
                            mod2a$coefficients['con.con_wamonth'], mod3a$coefficients['con.con_wamonth'], mod4a$coefficients['con.con_wamonth'], 
                            mod2a$coefficients['con.con_waday'], mod3a$coefficients['con.con_waday'], mod4a$coefficients['con.con_waday'],
                            0, 0, 0,
                            mod2a$coefficients['con.con_wbcouncil candidate'], mod3a$coefficients['con.con_wbcouncil candidate'], mod4a$coefficients['con.con_wbcouncil candidate'],
                            mod2a$coefficients['con.con_wborganizer'], mod3a$coefficients['con.con_wborganizer'], mod4a$coefficients['con.con_wborganizer'],
                            mod2a$coefficients['con.con_wbvoter'], mod3a$coefficients['con.con_wbvoter'], mod4a$coefficients['con.con_wbvoter'],
                            mod2a$coefficients['con.con_wbfriend'], mod3a$coefficients['con.con_wbfriend'], mod4a$coefficients['con.con_wbfriend'],
                            0, 0, 0,
                            mod2a$coefficients['con.con_wc_same'], mod3a$coefficients['con.con_wc_same'], mod4a$coefficients['con.con_wc_same'],
                            0, 0, 0,
                            mod2a$coefficients['con.con_wdbeaten'], mod3a$coefficients['con.con_wdbeaten'], mod4a$coefficients['con.con_wdbeaten'],
                            mod2a$coefficients['con.con_wdabducted'], mod3a$coefficients['con.con_wdabducted'], mod4a$coefficients['con.con_wdabducted'],
                            mod2a$coefficients['con.con_wdkilled'], mod3a$coefficients['con.con_wdkilled'], mod4a$coefficients['con.con_wdkilled'],
                            0, 0, 0,
                            mod2a$coefficients['con.con_wean opposition activist'], mod3a$coefficients['con.con_wean opposition activist'], mod4a$coefficients['con.con_wean opposition activist'],
                            mod2a$coefficients['con.con_weyour friend'], mod3a$coefficients['con.con_weyour friend'], mod4a$coefficients['con.con_weyour friend']
                ),
                'ses' = c(0, 0, 0,
                          mod2a$se['con.con_wamonth'], mod3a$se['con.con_wamonth'], mod4a$se['con.con_wamonth'], 
                          mod2a$se['con.con_waday'], mod3a$se['con.con_waday'], mod4a$se['con.con_waday'],
                          0, 0, 0,
                          mod2a$se['con.con_wbcouncil candidate'], mod3a$se['con.con_wbcouncil candidate'], mod4a$se['con.con_wbcouncil candidate'],
                          mod2a$se['con.con_wborganizer'], mod3a$se['con.con_wborganizer'], mod4a$se['con.con_wborganizer'],
                          mod2a$se['con.con_wbvoter'], mod3a$se['con.con_wbvoter'], mod4a$se['con.con_wbvoter'],
                          mod2a$se['con.con_wbfriend'], mod3a$se['con.con_wbfriend'], mod4a$se['con.con_wbfriend'],
                          0, 0, 0,
                          mod2a$se['con.con_wc_same'], mod3a$se['con.con_wc_same'], mod4a$se['con.con_wc_same'],
                          0, 0, 0,
                          mod2a$se['con.con_wdbeaten'], mod3a$se['con.con_wdbeaten'], mod4a$se['con.con_wdbeaten'],
                          mod2a$se['con.con_wdabducted'], mod3a$se['con.con_wdabducted'], mod4a$se['con.con_wdabducted'],
                          mod2a$se['con.con_wdkilled'], mod3a$se['con.con_wdkilled'], mod4a$se['con.con_wdkilled'],
                          0, 0, 0,
                          mod2a$se['con.con_wean opposition activist'], mod3a$se['con.con_wean opposition activist'], mod4a$se['con.con_wean opposition activist'],
                          mod2a$se['con.con_weyour friend'], mod3a$se['con.con_weyour friend'], mod4a$se['con.con_weyour friend']
                ))

p$lo <- p$diffs - 1.96*p$ses
p$hi <- p$diffs + 1.96*p$ses
p$loc <- rev(seq(1, dim(p)[1]))


pdf('03_Figures/scenchars_angerfearrally_coefs.pdf')
ggplot(p, aes(colour = Outcome)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, lwd = 2) +
  geom_point(aes(y = loc, x = diffs), size = 5) +
  geom_errorbarh(aes(y = loc, xmin = lo, xmax = hi),
                 lwd = 1) +
  scale_y_discrete(name ="Scenario Characteristic", limits = rev(seq(1,17)*3-1),
                   labels=unique(p$Characteristic)) +
  scale_x_continuous(name = "Coefficients") +
  theme_minimal() +
  theme(text = element_text(size = 12), 
        legend.position = 'bottom')
dev.off()




#####
## Table I.1: interaction with fear induction 
#####

con$treat_all <- car::recode(con$treat_assign, "'C'=0;else=1")

spec0 <- "con.con_wa_rec + con.con_wb_sim + con.con_wc_same + con.con_wd_sev + con.con_we_cred + se_pc_st + treat_all + community + surveyor_id"
spec1 <- " + party_close_st + female + age_st + education_st + assets_pc_urb_st + assets_pc_rur_st "
spec2 <- "+ activism_pc_st + trauma_pc_st"
spec3 <- " + con.con_wa_rec:treat_all + con.con_wb_sim:treat_all + con.con_wc_same:treat_all + con.con_wd_sev:treat_all + con.con_we_cred:treat_all"

## probability of rally

mod4a <- lm(as.formula(paste0('act_st ~ ', spec0)), data = con)
mod4a$se <- coeftest.cluster(con, mod4a, 'id')[,2]
mod4b <- lm(as.formula(paste0('act_st ~ ', spec0, spec3, spec1)), data = con)
mod4b$se <- coeftest.cluster(con, mod4b, 'id')[,2]
mod4c <- lm(as.formula(paste0('act_st ~ ', spec0, spec3, spec1, spec2)), data = con)
mod4c$se <- coeftest.cluster(con, mod4c, 'id')[,2]
mod4d <- lm(as.formula(paste0('act_st ~ ', spec0, spec1, spec2, spec3)), data = con)
mod4d$se <- coeftest.cluster(con, mod4d, 'id')[,2]


## anger - fear
mod2a <- lm(as.formula(paste0('ratio_af_st ~ ', spec0)), data = con)
mod2a$se <- coeftest.cluster(con, mod2a, 'id')[,2]
mod2b <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, spec3, spec1)), data = con)
mod2b$se <- coeftest.cluster(con, mod2b, 'id')[,2]
mod2c <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, spec3, spec1, spec2)), data = con)
mod2c$se <- coeftest.cluster(con, mod2c, 'id')[,2]
mod2d <- lm(as.formula(paste0('ratio_af_st ~ ', spec0, spec1, spec2, spec3)), data = con)
mod2d$se <- coeftest.cluster(con, mod2d, 'id')[,2]

stargazer(mod4a, mod4b, mod4c, mod2a, mod2b, mod2c, 
          se = list(mod4a$se, mod4b$se, mod4c$se, mod2a$se, mod2b$se, mod2c$se),
          covariate.labels = c('Time to Election', 'Non-Activist Victim', 'Same Province',
                               'Violence Severity', 'Source Credibility',
                               'Self Efficacy', 'Fear Treatment', 'Closeness to Party',
                               'Female', 'Age',
                               'Education', 'Urban Assets', 'Rural Assets',
                               'Past Activism', 'Past Trauma',
                               'Fear Treatment X Time to Election',
                               'Fear Treatment X Non-Activism Victim',
                               'Fear Treatment X Same Province',
                               'Fear Treatment X Violence Severity',
                               'Fear Treatment X Source Credibility',
                               'Constant'),
          omit = c('community', 'surveyor_id'),
          keep.stat = c('n', 'rsq'),
          no.space = T, digits = 2, star.cutoffs = c(.05, .01))


